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We introduce a randomized algorithm for computing the minimal-norm solution to 
an underdetermined system of linear equations. Given an arbitrary full-rank matrix 
Anxn with m < n, any vector b mx i, and any positive real number e less than 1, the 
procedure computes a vector x nx \ approximating to relative precision e or better the 
vector p n xi of minimal Euclidean norm satisfying A mxn Pnxi = b mx i- The algorithm 
typically requires 0(mn \og(y/n/e) + m 3 ) floating-point operations, generally less than 
the 0{m? n) required by the classical schemes based on Q-R-decompositions or bidiag- 
\ onalization. We present several numerical examples illustrating the performance of the 

' algorithm. 

m ' 

: 1 Introduction 

in 

Underdetermined systems of linear equations have arisen frequently in modern statistics 
and data analysis, and have been attracting much attention recently in various application 
domains; see, for example, [3], [I], and [6]. The solutions to underdetermined systems are 
not unique; the present article focuses on the solutions whose Euclidean norms are minimal. 

Given a full-rank matrix A mxn and a vector b mxl , with m < n, we would like to compute 
an accurate approximation to the vector p nx \ of minimal Euclidean norm satisfying 

Classical algorithms using Qi?-decompositions or bidiagonalization require 

^classical = 0(lll 2 7l) (2) 

floating-point operations in order to compute p n xi (see, for example, Chapter 5 in [S]). 

The present paper introduces a randomized algorithm that, given any positive real num- 
ber e less than 1, computes a vector x nx i approximating p nx i to relative precision e or better 
with respect to the Euclidean norm, that is, the algorithm produces a vector x nx i such that 

|| 2-nxl Pnxi || — £ ||Pnxl II, (3) 
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where || • || denotes the Euclidean norm. This algorithm typically requires 

C r and = 0(mn ]og(y/n/e) + m 3 ) (4) 

floating-point operations. When m is sufficiently large and n is much greater than m (that is, 
the system of linear equations is highly underdetermined) , then the cost in (J3J) is less than the 
cost in ([2]). Moreover, in the numerical experiments of Section [7J, the algorithm of the present 
article runs substantially faster than the standard methods based on Qi?-decompositions. 

The present paper describes an algorithm optimized for the case when the entries of 
A mxn and 6 mx i are complex valued. Needless to say, real- valued versions of our scheme 
are similar. The present article has the following structure: Section [2] sets the notation. 
Section [3] discusses a randomized linear transformation which can be applied rapidly to 
arbitrary vectors. Section H] describes the algorithm of the present paper. Section [5] proves 
that the procedure succeeds with high probability. Section [6] estimates the computational 
costs of the algorithm. Section [7J illustrates the performance of the scheme via several 
numerical examples. Section [8] contains several concluding comments. 



2 Notation 

In this section, we set notational conventions employed throughout the present paper. 

We abbreviate "independent and identically distributed" to "i.i.d." We consider the 
entries of all vectors and matrices in this paper to be complex valued. For any vector x, we 
define to be the Euclidean (/ 2 ) norm of x. For any matrix A, we define A* to be the 
adjoint of A. We define the condition number of A to be the I 2 condition number of A, that 
is, the greatest singular value of A divided by the least singular value of A. 

For any positive integer n, we define the discrete Fourier transform F nxn to be the matrix 
with the entries 

F jtk = -L e -^u-m-i)/n (5) 



for j, k = 1, 2, . . . , n — 1, n, where i = y/—l and e = exp(l) 



3 Preliminaries 

In this section, we discuss a subsampled randomized Fourier transform. [I], [7], [10], and [11] 
introduced a similar transform for similar purposes (these articles motivated us to write the 
present paper). 

For any positive integers / and n with / < n, we define the / x n SRFT to be the random 
matrix 

Tlxn = Glxn H nxn , (6) 

where Gi xn and H nxn are defined as follows. 

In dH]), Gi xn is the random matrix given by the formula 

G[xn = Si xn F nxn D nX ni (7) 
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where Si xn is the matrix whose entries are all zeros, aside from a single 1 in column Sj of 
row j for j = 1, 2, . . . , I — 1, /, where Si, s 2 , • • • , sz-i, s/ are i.i.d. integer random variables, 
each distributed uniformly over {1,2, ... ,n — l,n}; moreover, F nxn is the discrete Fourier 
transform defined in (jSJ), and D nxn is the diagonal matrix whose diagonal entries d\, d 2 , . . . , 
d n -i, d n are i.i.d. complex random variables, each distributed uniformly over the unit circle. 
(In our numerical implementations, we drew s±, s 2 , Sz_i, s/ from {1, 2, . . . , n — 1, n} 
without replacement, instead of using i.i.d. draws.) We observe that both F, 
are unitary. 



,, yn and D nxn 



In (ED 



, H nxn is the random matrix given by the formula 



H nxn @nxn rinxn ^nxn Onxn rinxri ^nxnj 



where II nxn and II nxn 
random, and Z nxn 



are permutation matrices chosen independently and uniformly at 
and 4xn are diagonal matrices whose diagonal entries Ci> C2> • • • > Cn-i, Cn 
and Ci, C2, ■ • • , Cn-i> Cn are i-i-cL- complex random variables, each distributed uniformly over 
the unit circle; furthermore, Q nxn and Q n xn are the matrices defined via the formulae 
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and (the same as but with tildes) 



/ cos(^) sin(£i) \ 
-sin(^i) cos(0 x ) 



a 



o 
o 





10 

'•• 
1 / 

/ 1 \ 

cos(# 2 ) sin(# 2 ) 

-sin(0 2 ) cos(# 2 ) 

1 

\ '•• j 

( '■■ o o o o N 

10 

cos(fl n _ 2 ) sin(fl n _ 2 ) 

-sin(0 n _ 2 ) cos(# n _ 2 ) 

\ 1 j 

( 1 

'•• 
1 

cos(£ n _i 







\ 






sin(6» n _i) 

V - sin(^ n _ 1 ) cos(0 n _i) j 



(10) 



where 9%, # 2 , ■ ■ ■ , # n _ 2 , 9 n -x, 9i, 9 2 , . . . , # n _ 2 , 9 n -\ are i.i.d. real random variables drawn 
uniformly from [0,2tt]. We observe that 6 nxn , Q nxn , n nxn , fl nxn , Z nxn , and Z nxn are all 
unitary, and so H nxn is also unitary 

We call the transform T) xn an "SRFT" for lack of a better term. 

The following technical lemma is a slight reformulation of Lemma 4.4 of 



Lemma 3.1. Suppose that a and 3 are real numbers greater than 1, and I, m, and n are 
positive integers, such that 

a 2 3 9 . 

n>l>- ^—m 2 . 11 

[a- \y 

Suppose further that Ti xn is the SRFT defined in (0|), and that Q nxm is a matrix whose 
columns are orthonormal. 

Then, the least singular value a m of Ti xn Q nxm satisfies 



o- m > \l — 




l_ 

an 



(12) 



with probability at least 1 



P m 



4 



4 Description of the algorithm 



Suppose that e is a positive real number less than 1, and I, m, and n are positive integers 
with m < I < n. Suppose further that A mxn is a full-rank matrix, b mx i is a vector, and p n xi 
is the vector of minimal Euclidean norm satisfying A mxn p nx \ — b mX \. In order to construct 
a vector x nx \ such that \\x nx i — p nx i\\ < £ lbnxi|| with high probability (increasingly high 
probability as a parameter a > 1 increases), we compute the vector c nx \ of minimal Euclidean 
norm that is a linear combination of random vectors and satisfies A mxn c nx \ = b mx i, then 
use the algorithm of [9 J to compute the orthogonal projection of c nx i onto the column span 
of (A mxn )*. More precisely, we perform the following five steps: 

1. Construct the matrix 

Slxm = Ti xn (A mxn ) , (13) 

applying the SRFT T\ xn defined in fl6]) to every column of (A mxn )* (see, for example, 
Subsection 3.3 of [12] for details on applying the SRFT rapidly). 

2. Construct the vector z\ x \ of minimal Euclidean norm solving the system of linear 
equations 

(Si xm) ^Ixl b m x i , (14) 

where Si xm is the matrix defined in f|T3|) (see, for example, Algorithm 5.7.2 in [8] for 
details on constructing Z[ X \). 

3. Construct the vector 

xn) z lxl, 

(15) 

where z% x \ is the vector of minimal Euclidean norm solving (|14p . and T\ xn is the same 
realization of the SRFT as in f|T3|) (see, for example, Subsection 3.3 of (12] for details 
on applying the adjoint of the SRFT rapidly). 

4. Use the algorithm of [9] for the construction of a vector y m xi minimizing 

mxn) Vm X 1 C n x 1 II 2 (16) 

to relative precision e 2 l/(an) or better, where c nx i is the vector defined in ( |T5l) . The 
parameter / for the algorithm of [9] should be the same as for the present algorithm. 

5. Construct the desired vector 

^nxl = (^mxn) Vmxlj (1^) 

where y m x\ is the vector from Step 4. 

Remark 4.1. In Step 2 above, we assume that Si xm defined in ( 1131) is a full-rank matrix. 
Lemma [3. II in Section [3] above guarantees this with high probability when / > m 2 , by taking 
a in the lemma arbitrarily large; numerical experiments indicate that / > m suffices. 

Remark 4.2. It is possible to improve the approximation x nx \ via preconditioned conjugate 
gradient iterations similar to those proposed in [9] . However, the approximation produced by 
the above algorithm is already highly accurate (see, for example, Section [7] or Theorem 15.41 
in Section O below), and further iterative improvement may double the running time of the 
algorithm. 
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5 Proof of accuracy 



In this section, we prove Theorem 15 .41 guaranteeing that the algorithm of Section H] produces 
high accuracy with high probability. 

The following lemma states that the orthogonal projection onto the column span of 
(Anxn)* of the vector c nx i defined in (TT51) is the vector of minimal Euclidean norm that the 
algorithm aims to approximate. 

Lemma 5.1. Suppose that I, m, and n are positive integers with m < I < n. Suppose further 
that A mxn is a full-rank matrix, 6 mx i is a vector, Si xm is the matrix defined in (T73|) and is 
a full-rank matrix, and c nx i is the vector defined in l[T5\) . 

Then, the orthogonal projection p nx i of c nx i onto the column span of (A mxn )* is the 
vector of minimal Euclidean norm satisfying 

A m xn Pnxl = ®mxl- (18) 

Proof. Combining ljI5j) . (IT41. and (fT5) yields that 

^mxnCnxl = ^mxl- (1^) 

Combining (Tl9j) and the fact that p nx \ is the orthogonal projection of c nx i onto the column 
span of (A mxn )* completes the proof. □ 

The following lemma states that, with high probability, the Euclidean norm of the vector 
c nx i defined in ffTol) is not too much greater than the Euclidean norm of its orthogonal 
projection p nx \ onto the column span of (A mxn )*. 

Lemma 5.2. Suppose that a and (3 are real numbers greater than 1, and I, m, and n are 
positive integers, such that (Tl\) holds. Suppose further that A 

mxn is a full-rank matrix, b mx \ 

is a vector, c nx i is the vector defined in ( T731) . and p nx \ is the orthogonal projection of c nx i 
onto the column span of (A mxn )* . 
Then, 




an 

\Cnxl\\ < \j— \\Pnxl\\ (20) 



with probability at least 1 — 4. 



Proof. Using the fact that A mxn is a full-rank matrix, we construct a Q-R-decomposition 

(A mxn ) QnxmRmxm (21) 

such that the columns of Q nxm ar e an orthonormal basis for the column span of (^4 mX n)*- 
We first show that the SRFT T Lxn used in (fl3j) and (jT5l) provides 

INxlH < \j — || (Qnxm)* (Tlxn)* ^xl|| (22) 

with probability at least 1 — -|, where Zi x i is the vector of minimal Euclidean norm solv- 
ing (fl4l) . We then express the left- and right-hand sides of (|22l) in terms of c nx i and p nx i, 
rather than zi x i, in order to obtain (1201) . 
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It follows from the fact that zixi is the vector of minimal Euclidean norm solving ([I] 
that zixi belongs to the column span of Si xm from (TH[) . that is, there exists a vector w mx i 
such that 

Zlxl = Si 

x m W m x l • (23) 

Combining ([23]), ([131) , and (EE} yields that 

H^Jxl IP = { z lxl)* Zlxl = (Wmxl)* (Rmxm)* {Qnxm)* {Tlxn)* z lxl- (24) 

The Cauchy-Schwarz inequality yields that 

mxm nxm xn) 2-ixl — -Rrnxm^mxl II II (Q nxm I Xn )*^xi||- (25) 
It follows from ([T2j) that 

1 1 Rmxm ^mxl || — W -J Wlxn Qnxm Rmxm W m xl\\ (^6) 

with probability at least 1 — 4. Combining f[2"T[) . ([TBI , and f[2"3"]) yields that 

Tlxn Qnxm Rmxm "'mxl = ^2x1- (27) 



Combining ([21]), $2>D, and (T27D yields 

We now express the left- and right-hand sides of ([22]) in terms of c nx i and p nx i, rather 
than zixi- 

First, we consider the left-hand side of ([22]) . Combining (fl5[) and the fact that the 
columns of (TJ Xn )* are orthonormal yields that 

II^XlH = ||Cnxl||- (28) 

Next, we consider the right-hand side of f[2~2"]) . It follows from the fact that the columns 
of Q n xm are an orthonormal basis for the column span of (A mxn )* that the orthogonal 
projection p nx i of c nx i onto the column span of {A mxn )* is 

Pnxl = Q nxm {Qnxm) Cnxi- 

(29) 

Combining ([29]) and the fact that the columns of Q n xm are orthonormal yields that 

||Pnxl|| = H(Qnxra) Cnxl||- 

(30) 

Combining ([30]) and ([T5]) yields that 

Ibnxlll = || {Qnxm)* {Tlxn)* Zlxl ||- (31) 

Finally, combining ([22]), ([25]), and ([51]) yields ([20]). □ 

The following lemma states that the vector x nx i produced by the algorithm is an accurate 
approximation to the orthogonal projection onto the column span of {A mX n)* of the vector 
c n xi defined in (FIB"]) , provided that the projection p nx i satisfies ( [211 . 
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Lemma 5.3. Suppose that e and a are positive real numbers with e < 1 < a, and I, m, and 
n are positive integers with m < I < n. Suppose further that A mxn is a full-rank matrix, 
bmxi is a vector, Si xm is the matrix defined in and is a full-rank matrix, c nx i is the 
vector defined in < f?3)) . p n xi is the orthogonal projection of c nx i onto the column span of 
(i mxn )*, Vmxi is a vector minimizing l[Tb)) to relative precision e 2 1 / (an) or better, and x nx \ 
is the vector defined in (T7[). Suppose in addition that [2D\l holds. 



Then, 

\\Xnxl -Pnxill < £ Ibnxill- (32) 

Proof. It follows from (TITj) that x nx \ belongs to the column span of (A mxn )*. Combining this 
fact and the fact that c nx i — p nx i is the orthogonal projection of c nx \ onto the orthogonal 
complement of the column span of (A mxn )* yields that c nx \— p nx \ is the orthogonal projection 
of 

c nxi — x nxi onto the orthogonal complement of the column span of (A mxn )* . Similarly, 
Pnxi — %nxi is the orthogonal projection of c nx i — x nx i onto the column span of (A mxn )*. 
We thus obtain the Pythagorean identity 

|| Cnxi Pnxi || ~t~ 1 1 Pnxi ^nxl || = || c nxl ^nxl || • (^3) 

It follows from the fact that p n xi is the orthogonal projection of c„ x i onto the column 
span of (A mxn )* that the minimal value of f|T6|) is \\p nx \ — c nx i|| 2 . Combining this fact, f|T7|) . 
and the fact that y mx i minimizes f lT6l) to relative precision e 2 l/(an) or better yields that 

2 2 ^ 2 

H^nxl c nxl|| ||Pnxl C nx i|| < ||Pnxl C nx i|| . (34) 

an 

Combining ( l33l) and ( |34l) yields that 

2 e 2 * 2 

H^nxl Pnxl|| ^ || c nxl Pnxl|| • (35) 

an 

It follows from the fact that c nx i — p n xi is the orthogonal projection of c nx i onto the 
orthogonal complement of the column span of (A mxn )* that 

||c n xl — Pnxl|| < ||c„xl||- (36) 

Combining (155] . (1561). and (BD|1 yields (J32]). □ 

Combining Lemmas 13. 1[ 15. 1[ 15.21 and 15.31 yields the following theorem, guaranteeing that 
the algorithm produces high accuracy with high probability. 

Theorem 5.4. Suppose that e, a, and (3 are positive real numbers with e < 1 < a and 
(3 > 2, and I, m, and n are positive integers, such that [Tl\) holds. Suppose further that 
A mxn is a full-rank matrix, b mx i is a vector, and p nx \ is the vector of minimal Euclidean 
norm satisfying 

AmxnPnxl = ^mxl- (37) 

Then, the vector x nx i defined in (T7\) satisfies 



\\Xnxl - Pnxi || < £ \\P nxl II (38) 

with probability at least 1 — 1. 
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Remark 5.5. The probability of failure in Theorem 15.41 is 2/(3, rather than just the proba- 
bility 1 /(3 of failure in Lemma I5.2[ in order to account for the possibility that the algorithm 
of [H] fails to produce a vector y mx i minimizing (TIT)]) to relative precision e 2 l/(an) or better, 
after using at most the number of floating-point operations specified in Section [6] below. 

Remark 5.6. Empirically, requiring (fTTj) appears to be excessive. In practice, choosing any 
I > 4m and a = 4 makes failure of the algorithm too improbable to detect (cf. Remark 2 
in [9j). 

6 Computational costs 

In this section, we tabulate the numbers of floating-point operations required by the five 
steps in the algorithm of Section HJ 

1. Applying TJ xn to every column of (A mxn )* costs 0{mn log(Z)). 

2. Constructing z\ x \ costs 0{lm 2 ). 

3. Applying (T ix „)* to z txl costs 0(n log(n)). 

4. Constructing ?// xl via the algorithm of [H] costs 0(mn log(an/e 2 ) + lm 2 ). 

5. Applying (A mxn )* to y mxl costs 0(mn). 

Summing up the costs in the five steps above, and using the facts that m < I < n and 
e < 1 < a, we see that the cost of the entire algorithm is 

^theoretical = 0(mn \og(^fan/e) + lm 2 ) (39) 

floating-point operations. Theorem 15.41 in Section [5] above guarantees that the algorithm 
produces high accuracy with high probability when / and a satisfy fllip . In practice, choosing 
I = 4m and a = 4 makes failure of the algorithm too improbable to detect (see also Remark 2 
in [9]), and the cost in (1391) then becomes 

Ctypkai = 0(mn log(vV^) + ™ 3 ) (40) 

floating-point operations. 



7 Numerical Results 

In this section, we describe the results of several numerical tests of the algorithm of the 
present paper. 

For various positive integers m and n with m < n, we use the algorithm to com- 
pute a vector x nx \ approximating the vector p nx \ of minimal Euclidean norm satisfying 
AmxnPnxi = &mxi, where b mx i is a vector, and A mxn is the matrix defined via the formula 

A m xn U mxm Smxm (Kxm) j (41) 



9 



in the experiments described below, U mxm is obtained by applying the Gram-Schmidt process 
to the columns of an m x m matrix whose entries are i.i.d. centered complex Gaussian random 
variables, V nxm is obtained by applying the Gram-Schmidt process to the columns ofannxm 
matrix whose entries are i.i.d. centered complex Gaussian random variables, and S mxm is a 
diagonal matrix, with the diagonal entries 

E jd = io-BU-V/(m-i) (42) 
for j — 1, 2, . . . , m — 1, m. Clearly, the condition number ka of A mxn is 

k a = E li i/E TO>m = 10 6 . (43) 
The vector b mx \ is defined via the formula 

bmxi = A mxn p nx i, (44) 
where p nx i is the vector defined via the formula 



Pnx 1 = —f= 6 3 V nl 1 > (45) 



m 

3=1 

with £]_, e 2 , . . . , e m _i, e m being pseudorandom positive and negative ones, and f£<i' t>ixi, 
• • • i v n^i l \ u ixi bein g tlie columns of V nxm . 

Remark 7.1. By construction, p nx i is the vector of minimal Euclidean norm such that 
AmxnPnxi = &mxij the Euclidean norm of p nx i is minimal since p nx \ belongs to the column 
span of (A mxn )*. The aim of the algorithm is to construct an approximation x nx \ to p n xi- 

For the direct computations, we used the classical algorithm for pivoted Q-R-decompo- 
sitions based on plane (Householder) reflections (see, for example, Chapter 5 in [S]). We 
implemented the algorithms in Fortran 77 in double-precision arithmetic, and used the La- 
hey/Fujitsu Express v6.2 compiler, with the optimization flag — o2 enabled. We used one 
core of a 1.86 GHz Intel Centrino Core Duo microprocessor with 2 MB of L2 cache and 
1 GB of RAM. For the algorithm of [9] used in Step 4 of the algorithm of Section @J we 
requested that y mX i minimize (TTBT) to relative precision (10~ 14 • ka) 2 4 m/n or better, where 
ka is the condition number of A mxn given in fT43|) . We used a double-precision version of 
P. N. Swarztrauber's FFTPACK library for fast Fourier transforms. 

Table 1 displays timing results with n = 16384 for various values of m; Table 2 displays 
the corresponding errors. Table 3 displays timing results with m = 256 for various values of 
n; Table 4 displays the corresponding errors. 

The headings of the tables have the following meanings: 

• m is the number of rows in the matrix A mxn , as well as the length of the vector & mx i, 
in A mxn p nx i = b mx \. 

• n is the number of columns in the matrix A mxn , as well as the length of the vector 
Pnxit in A mxn p nx i = b mx \. 
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Table 1 Table 2 



m 


n 


I 


to 




to/t T 


m 


n 


/ 


£o 




128 


16384 


512 


.27E1 


.24E1 


1.2 


128 


163t 


U 


512 


.14E- 


-14 


.16E- 


-14 


256 


16384 


1024 


.11E2 


.56E1 


2.0 


256 


163J 


J4 


1024 


.11E- 


-14 


.17E- 


-14 


512 


16384 


2048 


.60E2 


.20E2 


3.0 


512 


163J 


i4 


2048 


.80E- 


-15 


.29E- 


-14 



Table 3 Table 4 



m 


n 


I 


to 


tr 


to/t v 


m 


n 


I 


£o 




256 


4096 


1024 


.26E1 


.20E1 


1.3 


256 


4096 


1024 


.27E- 


-15 


.31E- 


-14 


256 


8192 


1024 


.51E1 


.32E1 


1.6 


256 


8192 


1024 


.45E- 


-15 


.27E- 


-14 


256 


16384 


1024 


.11E2 


.56E1 


2.0 


256 


16384 


1024 


.11E- 


-14 


.17E- 


-14 


256 


32768 


1024 


.29E2 


.16E2 


2.5 


256 


32768 


1024 


.22E- 


-14 


.16E- 


-14 



• / is the number of rows in the matrix TJ xn used in Steps 1 and 3 of the algorithm of 
Section HJ as well as the analogous parameter used in the algorithm of [9] needed in 
Step 4. 

• to is the time in seconds required by the direct, classical algorithm. 

• t T is the time in seconds required by the algorithm of the present paper. 

• t /t r is the factor by which the algorithm of the present paper is faster than the classical 
algorithm. 

• Eq is defined via the formula 

II (°) _ II 

Praxl Pnxl\\ 

£ ° = ii ii — ' 

K A ||Pnxl || 

where ka is the condition number of A mxn given in (j4"3"l) . and x^ ^ is the vector produced 
by the direct, classical algorithm approximating the vector p n xi of minimal Euclidean 
norm such that A mxn p nxl = b mxl . 

• e r is defined via the formula 

e r = \\^1_ZP^A } (47) 

K A ||Pnxl || 

where ka is the condition number of A mxn given in (j4"3"I) . and x nx \ is the vector produced 
by the algorithm of Section H] approximating the vector p nx \ of minimal Euclidean norm 
such that A mxn p nx i = b mx \. 

Remark 7.2. Standard perturbation theory shows that Eq and e r are the appropriately 
normalized measures of the relative precision produced by the algorithms; see, for example, 
Section 5.5.3 in [5]. 
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The values for e r reported in the tables are the worst (maximal) values encountered 
during 10 independent randomized trials of the algorithm, as applied to the same matrix 
A mxn and vector b mxl . The values for t T reported in the tables are the average values 
over 10 independent randomized trials. None of the quantities reported in the tables varied 
significantly over repeated randomized trials. 

The following observations can be made from the examples reported here, and from our 
more extensive experiments: 

1. When m = 512, n = 16384, and the condition number of A mxn is 10 6 , the randomized 
algorithm runs 3 times faster than the classical algorithm based on plane (Householder) 
reflections, even at full double precision. 

2. Our choice / = 4m appears to make failure of the algorithm too improbable to detect. 

3. The algorithm of the present article reliably produces high precision at reasonably low 
cost. 

8 Conclusion 

This article provides a fast algorithm for computing the minimal-norm solution to an under- 
determined system of linear equations. If the matrices A mxn and (A mxn )* associated with 
the system of linear equations can be applied sufficiently rapidly to arbitrary vectors, then 
the algorithm of the present paper can be accelerated further. 

The theoretical bounds in Lemma [37T1 Lemma |5T2| and Theorem 15 . 41 should be considered 
preliminary. Our numerical experiments indicate that the algorithm of the present article 
performs better than our estimates guarantee. Furthermore, there is nothing magical about 
the subsampled randomized Fourier transform defined in ([6]). In our experience, several other 
similar transforms seem to work at least as well, and we are investigating these alternatives 
(see, for example, [2]). 
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